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Statistics of Poincare recurrences is studied for the base-pair breathing dynamics of an all-atom 
DNA molecule in realistic aqueous environment with thousands of degrees of freedom. It is found 
that at least over hve decades in time the decay of recurrences is described by an algebraic law 
with the Poincare exponent close to /3 = 1.2. This value is directly related to the correlation decay 
exponent z/ = — 1, which is close to v ~ 0.15 observed in the time resolved Stokes shift exper¬ 

iments. By applying the virial theorem we analyse the chaotic dynamics in polynomial potentials 
and demonstrate analytically that exponent (3 — 1.2 is obtained assuming the dominance of dipole- 
dipole interactions in the relevant DNA dynamics. Molecular dynamics simulations also reveal the 
presence of strong low frequency noise with the exponent r] = 1.6. We trace parallels with the 
chaotic dynamics of symplectic maps with a few degrees of freedom characterized by the Poincare 
exponent [3 ~ 1.5. 
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The celebrated Poincare recurrence theorem of 1890 
[T] guarantees that a dynamical trajectory with a fixed 
energy and bounded phase space will always return in a 
close vicinity of the initial state. For dynamical systems 
with hard chaos the statistics of Poincare recurrences, 
and the related probability to stay in a bounded phase 
space region behave similarly to coin flipping and drop 
exponentially with the return time r mm- However, in 
the generic case of chaos with divided phase space, when 
islands of integrable motion are embedded in a chaotic sea 
HIS], it was established that the probability distribution 
of recurrences P{r) is described by an algebraic decay 

P(t) cx l/r'^ , (1) 

with the Poincare exponent /3 ^ 1.5 ISHII]. This slow 
decay originates from sticking of dynamical trajecto¬ 
ries in the vicinity of stability islands and results in a 
slow decay of the corresponding atuocorrelation function 
C{t) oc tP(t) [To]. Most of the studies of the Algebraic 
Statistics of Poincare Recurrences (ASPR) considered 2D 
symplectic maps, notably, the Chirikov standard map, 
with recurrence times changing by more than 10 orders 
of magnitude m- A few recent studies of Hamiltonian 
systems with a larger number of degrees of freedom also 
revealed an algebraic decay of recurrences with similar 
values of the Poincare exponent (3 ^ 1.3 — 1.5 [T3UT6] . 

The generic nature of the ASPR phenomenon is well 
established. It is known to occur on a huge range of phys¬ 
ical scales from electron trajectories for microwave ioniza¬ 
tion of Rydberg atoms [num to comet orbits in the Solar 
System m- One should expect that it is also inherent 
in conformational dynamics of macromolecules. These 
systems are characterized by complex energy landscapes. 


with numerous barriers and saddle points crossed during 
thermal motion. It is tacitly assumed that such dynam¬ 
ics results in a developed chaos with multi-exponential 
relaxation decay. Recent experimental evidences suggest, 
however, that the ASPR may play an important role, no¬ 
tably, in behavior of the double helical DNA. The power- 
law relaxation in the B-DNA double helix was discov¬ 
ered [20l [21] and carefully studied during the last decade 
by using the Time Resolved Stokes Shift (TRSS) exper¬ 
iments [22H22]. In the last years these results were an¬ 
alyzed with different theoretical approaches [28f[34] . but 
nevertheless, the possible underlying molecular mecha¬ 
nism remains elusive. By analogy with the long known 
power-law kinetics in proteins [35], this effect in DNA 
is interpreted in terms of models developed earlier for 
glassy systems, with multiple substates, hierarchical re¬ 
laxation, mode coupling, etc. However, unlike proteins, 
the B-DNA molecule has only a few well-studied confor¬ 
mational substates with relatively fast and spatially lo¬ 
calized dynamics that does not resemble those in glasses. 
With hydration water and surrounding ions included, the 
system becomes more complex, but, in spite of all efforts, 
there is no agreement even on whether the power-law 
relaxation in the sub-microsecond time range is due to 
DNA itself or the hydration water, or both [27l|32j|36]. 

We asked, what if, instead of the spin-glass like effect 
of multiple degrees of freedom, the power-law relaxation 
in DNA represents a manifestation of the ASPR phe¬ 
nomenon. According to experiment, the decay of the 
TRSS signal in DNA is described by time autocorrela¬ 
tion function C{r) oc 1/r^ with u ^ 0.15 [21], that is, 
rather close to ASPR in systems with a few degrees of 
freedom. This decay is observed over six decades in time 
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from 10“^^ to 10“^ sec, therefore, one can hope to detect 
it by means of all-atom Molecular Dynamics (MD) sim¬ 
ulations. The TRSS signal is obtained by substituting a 
polarity sensitive dye coumarine for one of the stacked 
bases, and the measured effect can be due to any motion 
that changes the dye local electric field 120 ]. Since in 
MD different motions are coupled we started from a trial 
search of relevant parameters. The statistics of Poincare 
recurrences for kinetic energies of selected atoms revealed 
only exponential decay. A similar behavior, except for 
minor details, was found for conformational transitions 
of backbone torsion angles. After a few unsuccessful tri¬ 
als, however, the ASPR has been revealed in the base-pair 
breathing motion. 




FIG. 1: (Color online) The left panel shows the Watson-Crick 
base pair formed by guanine and cytosine. The three H-bonds 
are shown by thick dashed lines. Their lengths are about 3 A. 
The right panel shows a snapshot of the model system taken at 
the end of one of the 65 MD trajectories involved in statistical 
analysis. A tetramer fragment of a B-DNA double helix was 
built from two identical strands with the self-complementary 
base pair sequence GCGC. The six sodium ions necessary for 
charge neutralization are shown as spheres. It is seen that the 
four-level stack of base pairs shown in the left panel remains 
stable during simulations. 

The base-pair breathing occurs due to temporal break¬ 
ing of one or more hydrogen bonds (H-bonds) in a 
Watson-Crick (WC) pair. The statistics of Poincare re¬ 
currences was studied for the model system shown in Fig. 

On the left, a GC pair is displayed with three H- 
bonds. The right panel shows a tetramer duplex formed 
by two strands with identical GC-alternating sequences. 
This is a minimal symmetrical structure with two exter¬ 
nal and two internal base pairs. The duplex is placed in 
a small water box of 489 water molecules with periodic 
boundaries. Six sodium ions are added for neutralization. 
All-atom MD simulations were carried out in internal co¬ 
ordinates, with fixed backbone bond lengths and rigid 
bases, using Hamiltonian equations and a symplectic 
integrator [38] with the time step of 0.01 ps [39]. The re¬ 
cent version of the AMBER force field [4Qlj4^ was used 
with SPC/E water [43]. The system had 3226 degrees of 


freedom for 1705 atoms. 

The base-pair breathing was followed by measuring the 
distances (R) between the H-bond forming atoms shown 
in Eig. [^at every time step. The stopwatch was started 
when a given distance exceeded a certain threshold {Rth) 
and stopped once the boundary was crossed in the op¬ 
posite direction. These events are called Poincare recur¬ 
rences. The integrated probability distribution P{r) is 
obtained by counting the number of recurrences with du¬ 
ration larger than r and normalizing it by the total num¬ 
ber of events, that is, P(0) = 1 by construction. Eunction 
P(r) is a very powerful instrument of analysis because it 
is positive definite and, due to statistical averaging over 
a large number of crossings, stable with respect to fluc¬ 
tuations (see e.g. discussion in [lO]). Importantly, these 
computations are trivially parallelizable, that is, the P{r) 
statistics can be accumulated in a large number of inde¬ 
pendent MD trajectories. Most of the results discussed 
below were obtained by using parallel computations on 
65 cores. 



FIG. 2: (Golor online) The left panel shows statistics of 
Poincare recurrences P(r) for the three Watson-Grick H- 
bonds shown in Fig. [^ The results for distances 06N4, 
N1N3, and N202 are displayed by solid black, dashed red 
and dotted blue curves, respectively. The threshold distance 
is Rth = 3.15 A in all three cases. The right panel shows the 
corresponding dependencies of average bond distances {R) ob¬ 
tained for recurrences of different duration. Here and in other 
figures the logarithms are decimal. 

Representative results obtained for the three H-bonds 
in terminal base pairs are shown in Fig. [^ The left panel 
displays double-logarithmic plots of P{r) distributions 
obtained with Pt/i=3.15 A. This threshold is close to the 
equilibrium H-bond lengths, therefore, a large fraction 
of recurrences result from oscillations within the bonded 
ground state. These motions give for r<0.3 ps a charac¬ 
teristic fall of P{r) typical for exponential decays. The 
right panel of Fig. [^ shows the plots of average distances 
for recurrences of different duration. For returns shorter 
than 0.3 ps these values remain around 3.5 A, that is, 
the H-bonds are not broken. With t>0.3 ps the aver¬ 
age distances grow with r and the P{r) decay becomes 
algebraic. The short time boundary of the power-law re¬ 
laxation is close to that in TRSS experiments [21]. More¬ 
over, the decay exponent A ^ 1 for all three H-bonds is 
not far from experimental /3 — 1 = u ^ 0.1b. To refine 
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these results we decided to concentrate upon the 06N4 
H-bond in terminal base pairs. This H-bond is broken 
easier than other, therefore, it is an adequate indicator 
of partial base-pair openings. The internal base pairs 
were also checked, but they are opened rarely and only 
traces of algebraic decay were detectable in the tails of 
P(r) distributions. 



FIG. 3: (Color online) The left panel shows statistics of 
Poincare recurrences P(r) for H-bond 06N4 measured at four 
different threshold distances Rth- The results for Pt/i=3.15, 
3.55, 4.15, and 4.55 A are shown by the solid black, dashed 
red, dotted blue, and thin solid magenta curves, respectively. 
The right panel shows similar statistics for the threshold of 
3.55 A evaluated in two separate runs. The solid black curve 
corresponds to the original simulation (dashed red curve in 
the left panel). The second run, shown by dashed red curve, 
is carried out with a flat-bottom restraint that prevented the 
06N4 distance to go below 3.15 A thus pushing the system to 
better sampling of long returns. In both panels the dash-and- 
dot straight line shows the power law decay with the Poincare 
exponent /3 = 1.15. 

To get a better estimate of the short-time boundary 
of the algebraic decay in P{r) we need to remove the 
contribution of quick returns that do not result in H-bond 
opening. To this end the threshold Rth was gradually 
increased up to 4.55 A. The results of these computations 
are displayed in the left panel of Fig. It is seen that 
the short time hump is essentially removed already with 
Pt/i=3.55 and that for the true base-pair breathing the 
algebraic decay starts from very short recurrences of only 
0.1 ps. This refined boundary is in excellent agreement 
with the experimental data m- 

These computations were also used to refine the es¬ 
timate of the exponent (3. Linear fits of the plots in 
the left panel of Fig. were carried out in the range 
—0.5 < log{r/ps) < 2.5, which gives exponent values 
/d=1.27, 1.20, 1.13, and 1.05 for thresholds i?t/i=3.15, 
3.55, 4.15, and 4.55 A, respectively. It is seen that there 
is a moderate influence of Rth upon the exponent, which 
can be related to the fact that longer recurrences corre¬ 
spond to larger atom-atom separations (Fig. [^. Note 
that the p values obtained give the exponent of correla¬ 
tion decay v = ^ — 1 very close to the experimental TRSS 
value (z/ ~ 0.15) [21]. For visual comparison, the P(r) 
decay predicted from experiment is shown in Fig. iby 
the dash-and-dot lines. 

The experimental long-time boundary occurs at «100 


ns and it corresponds to the maximal resolution of the 
TRSS method [50l [21]. The long-time boundary in our 
computations is limited by sampling. We tried to push it 
somewhat further by taking into account that H-bonds 
in our simulations spend almost all the time oscillating 
around the ground energy minimum. To reduce this non¬ 
productive time a fiat-bottom restraint was added that 
prevented the 06N4 distances to go below 3.15 A. It is 
understood that this simple ad hoc trick perturbs realis¬ 
tic dynamics, but the information it provides may be use¬ 
ful. The results of these computations shown in the right 
panel of Fig. [^confirm that the sampling is indeed im¬ 
proved, with the time range of the approximate power law 
extended by about an order of magnitude. For shorter 
duration the P{t) distribution reproduces the curve ob¬ 
tained without the restraint. Therefore, one can reason¬ 
ably expect that the long-time boundary of the power-law 
decay of P{t) occurs at least at ^10 ns. 
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FIG. 4: (Golor online) The average power spectral density 
of time fluctuations of the 06N4 distance obtained from 16 
independent trajectories of 10 ns each, with the data stored 
at every time step. The dash-and-dot straight line shows the 
power law decay with the exponent rj = 1.63 obtained from 
the linear fit in the range — 5 < log(/ • ps) < —1. The in¬ 
sert shows a representative time trace of the same parameter 
between two consecutive crossings of the threshold distance 
Rth = 3.15 A. 

Encouraged by the foregoing findings, notably, the sur¬ 
prisingly good agreement with experiment obtained with 
no adjustable parameters, we decided to study the power 
spectrum of the fluctuations of the 06N4 distance. The 
results shown in Fig. [^reveal a strong power-law growth 
of density S for low frequencies. However, the exponent 
T] = 1.63 estimated by linear fitting does not seem to be 
related to neither experimental data nor the Poincare re¬ 
currences studied above. This paradoxical observation is 
discussed further below. 

To shed light on the possible origin of the foregoing re¬ 
sults we use the virial theorem m for analysis of chaotic 
dynamics in systems with representative polynomial po- 
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tentials. The corresponding Hamiltonian reads 

H{p^r) = /2 — a/r'^ = E ^ (2) 

where (p, r) is the pair of conjugated variables, with 
momentum p and radial coordinate r. The energy 
^ < 0 corresponds to a bounded motion; a ^ 1 and 
m > 0 are numerical coefficients (the mass is taken as 
unity). We assume that in addition to radial dynam¬ 
ics there is a chaotic motion in angle degrees of free¬ 
dom. For this system the action can be estimated as 
J pr (5^)(^-2)/(2m) ^ certain numerical con¬ 

stant 6 < 0, which follows from relationships p^ ^ 
and E ^ p^ . Thus we have bE ^ \E\ ^ j 2 ml{m- 2 )^ 
uj = dH/dJ ^ j(m+2)/(m-2) ^ where cj is a fre¬ 
quency of motion and r is the related characteristic time 
scale. 

As a result, the measure p related to the sticking 
time scale r is obtained as /i ^ J ^ i/^(’^- 2 )/(m+ 2 )^ 
which follows from the fact that in Hamiltonian systems 
the measure is proportional to the phase volume, that 
is, /i ^ J JdO ^ J X 27 t^ where 0 is the angle vari¬ 
able conjugated to action J. According to HO] we have 
P(r) ^ dp/dr ^ p/r where /r(r) is the measure of a 
region where a trajectory is stuck for the time r. This 
follows from the ergodicity relation according to which 
the measure of a region is proportional to the time spent 
by the trajectory in this region pipr) ^ tP{t)I (r), where 
(r) = P{r)dr is the average time of recurrences OISj. 
From these relations one obtains the following expression 
for the Poincare exponent 

P = 2m/{m + 2). (3) 

Here P(r) can be considered as an integrated probability 
of Poincare recurrences or as a survival probability in a 
given region for time periods longer than r since both are 
proportional to each other [IQIIII]. 

Consider some earlier studied potentials with different 
m. With m = 1 we get the Kepler problem. This case 
appears in the microwave ionization of Rydberg atoms 
iniiiH], and also in the comet m or dark matter [45] 
dynamics in the Solar System affected by Jupiter. In 
both cases the energy change occurs when the particle 
passes near the perihelion (near the nuclei or near the 
Sun). This energy change produces chaotic dynamics in 
the system. In this case the measure p ^ J 1/ is 
diverging at \E\ ^0 and from (§ we have /3 = 2/3 < 1, 
which agrees with the analytical and numerical results re¬ 
ported in [46]. For m = 2 we have a period independent 
of action, which gives P = 1 without decay of correla¬ 
tions, that is, C ^ tP{t) oc const with = /3 — 1 = 0. 

For the most relevant case of dipole-dipole interactions 
that are dominant in H-bonds and, more generally, in 
neutral polar systems like B-DNA with ions in water, 
we have m = 3 and (3 = 1.2. The last value is close 
to that obtained in our numerical simulations as well as 


that corresponding to exponent n = ^ — 1 found in the 
TRSS experiments. We believe that the above estimates 
correctly capture the main physical effects in the dynam¬ 
ics of this complex system and are at the origin of the 
observed slow algebraic decay of Poincare recurrences. It 
is understood that in systems with thousands of atoms 
like that studied here other factors can interfere. It is 
possible that Coulomb forces (m = 1) of locally uncom¬ 
pensated charges are responsible for a certain reduction 
of /3 to a slightly smaller value compared to the above 
theoretical estimate. We also note that for the van der 
Waals potential we have m = 6 with the corresponding 
p = 3/2. 

Finally, consider the result shown in Fig. [^ The mea¬ 
sure of sticking regions is p r\j J r\j \E\^/ ^ r\j r\j 

1/r^/^. It decreases with large r, but the typical atom- 
atom separations are growing as r . The 

correlation function C(r) = (r(t + T)r{t))t is commonly 
defined for a bounded variable r(t), so that the average 
square variation is {r{r)^)t ^ t f C{r)dT ^ where 

u = P — 1 < lis the correlation decay exponent related to 
the super-diffusive growth. In the present case, we have 
a variable that grows with r, which gives an additional 
contribution to the average square variation (r(r)^) ^ 
^ ^ r'^ with k = 3p — 1 = 2.6. Using 

the Wiener-Khintchine relation (see e.g. 03 ) between 
the square variation of a time dependent variable r{t) 
and its spectral density S{f) we obtain 

S{f) oc 1/r , 7? = K - 1 = 3/3 - 2 = 1.6 (4) 

in good agreement with Fig. In other words, we have 
a very strong low frequency noise, with the exponent p 
larger than usual [48] , that results from only the inherent 
chaotic dynamics of the system, with no external noise 
involved. 

In summary, using all-atom MD simulations, we un¬ 
covered the existence of algebraic decay of Poincare re¬ 
currences in the base-pair breathing motion of the B- 
DNA molecule, with the decay exponent P ^ 1.2 and a 
strong divergence of spectral density of vibration motion 
with the exponent ^ 1.6 over at least five decades in 
time from 10“^^ to 10“^ sec. These results are well de¬ 
scribed by the proposed theory of Poincare recurrences 
for chaotic dynamics in polynomial potentials, assum¬ 
ing a dominant contribution of dipole-dipole interactions. 
The theory correctly captures the qualitative origin of 
this effect in MD in spite of the difference in the number 
of degrees of freedom. The exponent P ^ 1.2 predicts 
a power-law relaxation of the corresponding correlations 
with the exponent u = P — 1 ^ 0.2. This prediction as 
well as the time scale of the algebraic decay in MD are 
in striking similarity with TRSS experiments, suggesting 
that this effect is due to the base-breathing motion. Fur¬ 
ther studies should help to clarify the relationship be¬ 
tween the inherent dynamical chaos in DNA and these 
experimental data. 
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